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Restricted path integral Monte Carlo simulations are used to calculate the equilibrium properties 
of hydrogen in the density and temperature range of 9.83 x 10 -4 < p < 0.153 gem -3 and 5000 < 
T < 250 000 K. We test the accuracy of the pair density matrix and analyze the dependence on the 
system size, on the time step of the path integral and on the type of nodal surface. We calculate 
the equation of state and compare with other models for hydrogen valid in this regime. Further, 
we characterize the state of hydrogen and describe the changes from a plasma to an atomic and 
molecular liquid by analyzing the pair correlation functions and estimating the number of atoms 
and molecules present. 



I. INTRODUCTION 

In spite of the simple composition, hydrogen exhibits a 
surprisingly complex phase diagram, which is the subject 
of numerous experimental and theoretical approaches. 
In this work, we study the high temperature regime of 
5000 < T < 250 000 K where hydrogen undergoes a 
smooth transition with increasing temperature from a 
molecular fluid through an atomic regime and finally to 
a two component plasma of electrons and protons (see 
Fig. [l]). The properties of hydrogen in this regime are 
crucial for the evolution of stars and the characteristics 
of the Jovian planets. 

A variety of simulation techniques and analytical mod- 
els have been developed to describe hydrogen at low den- 
sity. This regime has been studied with chemical models 
jj], ||, [| that describe hydrogen as a mixture of inter- 
acting molecules, atoms, free protons and electrons. The 
chemical composition is determined by minimizing an ap- 
proximate free energy function constructed from known 
theoretical limits. In this paper, we focus on low and 
intermediate densities 9.83 x 10~ 4 < p < 0.153 gem -3 
corresponding to 14 > 7' s > 2.6, where one expects the 
chemical models to work well although the properties 
of hydrogen are determined by the complex interplay of 
long-range Coulomb forces leading to strong coupling and 
bound states as well as degeneracy effects. 

All these effects can also be described from first prin- 
ciples simulation. There are ab initio methods such as 
restricted path integral Monte Carlo simulations (PIMC) 
J|, ^|, ^) and density functional theory molecular dynam- 
ics (DFT-MD) f], gj. The focus of the work is to test the 
equation of state (EOS) derived from chemical models 
and the actual density-temperature limits of the validity 
of the chemical picture. Additionally, we provide data 
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FIG. 1: Phase diagram of hydrogen as a function of tempera- 
ture and density for the different regimes: plasma, the atomic 
and the molecular regime. The dash-dot lines indicate the 
approximate boundaries. The crosses indicate the parame- 
ters for which PIMC simulations have been performed, the 
solid lines are isobars, and the dashed lines represent contour 
lines of constant permutation probability of the electrons (as 
indicated on the line.) 



to determine the parameters of the free energy models. 
Chemical models are expected to become inaccurate in 
regions of high density where the lifetime of molecules 
reduces to a few molecular vibrations J8|. 

We present results from new, more accurate, PIMC 
simulations. First, we analyze the accuracy of the pair 
density matrix and the error of the "time" discretization. 
Then we analyze the finite size dependence of the derived 
EOS and discuss the fixed node errors by comparing re- 
sults from simulations using free particle (FP) nodes and 
variational nodes. Furthermore, we calculate pair corre- 
lation functions which we use in conjunction with a clus- 
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ter analysis to characterize the state of hydrogen at dif- 
ferent temperatures and densities. We use atomic units 
(lengths in Bohr radii and energies in Hartrees) through- 
out this work except where indicated otherwise. 



II. PATH INTEGRAL MONTE CARLO 
METHOD 

A. Restricted Path Integral 

The density matrix (DM) of a fermion system at tem- 
perature kgT = 1/(3 can be written as an integral over 
all paths Rt, 



/5 (R ,R /3 ;/?) = ^£(-l) p I dR t e- s ^l (1) 

Rq— *7- > R. ( 3 

R t stands for the entire paths of N particles in 3 di- 
mensional space R t = (ri t , . . . , rjvt) beginning at Ro and 
ending at VRp. V labels the permutation of the particles 
and (— l) v to its signature. For non-relativistic particles 
interacting with a potential V(R), the action of the path 
S[Rt] is given by, 



S[R t 



dt 



dR{t) 



hdt 



V(R(t)) 



const. 



(2) 



In practice one discretizes H the path into a finite num- 
ber of imaginary time slices M corresponding to a time 
step t = /3/M. 

For fermionic systems the integration is complicated 
due to the cancellation of positive and negative contribu- 
tions to the integral, (the fermion sign problem). It has 
been shown (l(], [llj] that one can evaluate the path inte- 
gral by restricting the path to only specific positive con- 
tributions. One introduces a reference point R* on the 
path that specifies the nodes of the DM, p(R, R*, t) = 0. 
A node-avoiding path for < t < (3 neither touches nor 
crosses a node: p(R(t),R* ,t) ^ 0. By restricting the 
integral to node-avoiding paths, 



p F (Rp,R*;{3) = 

JdR Pf (Ro,R*;0) I dR t e ~ s[Rtl 

Ro^RaeT(R') 



(3) 



(T(R*) denotes the restriction) the contributions are 
positive and therefore PIMC represents, in principle, a 
solution to the sign problem. The method is exact if 
the exact fermionic DM is used for the restriction. How- 
ever, the exact DM is only known in a few cases. Most 
applications have approximated the fermionic DM by a 
determinant of single particle DMs, 



p(R,R';/3) = 



pi(ri,r' i; /3) 
pi(ri,r^;/3) 



pi(r w ,ri;$) 
pi(rAr,r^;/3) 



(4) 



This approach has been extensively applied using the 
free particle (FP) nodes derived from the single-particle 
density matrix 

pi(r, r',/3) = (4^ A/3)- 3 / 2 exp {-(r - r') 2 /4A/3} (5) 

with A = Ti 2 /2m, including applications to dense hydro- 
gen 0, |[ U . It can be shown that for temperatures larger 
than the Fermi energy, the interacting nodal surface ap- 
proaches the FP nodal surface. In addition, in the limit 
of low density, exchange effects are negligible: the nodal 
constraint has a small effect on the path and therefore 
its precise shape is not important. The FP nodes also 
become exact in the limit of high density when kinetic 
effects dominate over the interaction potential. However, 
for the densities and temperatures under consideration, 
interactions could have a significant effect on the nodal 
surfaces. 

To gain some quantitative estimate of the possible ef- 
fect of the nodal restriction on the thermodynamic prop- 
erties, it is necessary to try an alternative. In addition 
to FP nodes, we used nodal surface of a variational den- 
sity matrix (VDM) |i~2| derived from a variational princi- 
ple that includes interactions and atomic and molecular 
bound states. We assume a trial DM with parameters qi 
that depend on imaginary time [3 and R', 



p(R,R';/3) =p(R,gi, 



By minimizing the integral 
fdp(R,R';(3) 



dR 



di3 



Hp(R,R';j3)) =0 



(G) 



(7) 



one determines equations for the dynamics of the param- 
eters in imaginary time: 



1 dH <-> • 

- -rr^ + M q = where H 

2 aq 



pHp dR . (8) 



The norm matrix is: 

d 2 



J\fij = lim ■ 

q'^q , 



J dR p(R, q ; (3) p(R,q';(3) 



(9) 



We assume the DM is a Slater determinant of single par- 
ticle Gaussian functions 

pi(r, r',/3) = (7rw)- 3/2 exp{-(r-m)7w + d} (10) 

where the variational parameters are the mean m, 
squared width w and amplitude d. The differential equa- 
tions for this ansatz are given in |l2[ . The initial condi- 
tions at (3 — > are w = 2/3, m = r' and d = in order 
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to regain the correct FP limit. It follows from Eq. |t] 
that at low temperature, the VDM goes to the lowest en- 
ergy wave function within the variational basis. For an 
isolated atom or molecule this will be a bound state, in 
contrast to the delocalized state of the FP DM. A further 
discussion of the VDM properties is given in 0]. Note 
that these nodes are only used to determine the nodal 
restriction in Eq. (II. 3). The complete potential is taken 
into account in the path integral action as discussed in 
detail in 101. 



B. Accuracy of the Method 

The numerical implementation of the PIMC method 
requires one to make several approximations. Inaccura- 
cies can be caused by statistical errors from the MC inte- 
gration, inaccuracies in the numerically determined pair 
density matrices, a dependence on the time step of the 
path integral because of ./V-body (N > 3) correlations, 
finite size effects, and nodal errors from approximations 
in the trial density matrix. Their effects on the accuracy 
of the computed thermodynamic averages are quantita- 
tively estimated in this section. 

Statistical errors in the estimators for the thermody- 
namic quantities are calculated from the block averages 
generated by the MC simulations. The correlations be- 
tween blocks are taken into account by performing a 
blocking analysis |l3| . The resulting error bars (one stan- 
dard deviation) are given for all observables throughout 
this paper as a number in parenthesis referring to the 
least significant digit. 

In the following discussion, we compare the internal 
energy and the pressure calculated using the virial theo- 
rem: 



3Pv = 2K + V 



(11) 



where v is the volume of the simulation cell, K the kinetic 
energy and V the potential energy. Accurate estimation 
of the pressure requires a high accuracy in the kinetic 
and the potential energies because they tend to cancel 
out. In the molecular regime at low density, both terms, 
dominated by intra-molecular contributions, cancel to a 
large extent leaving behind the molecular gas pressure, 
which is much less than the inter-molecular forces. As 
a result, the pressure is, in general, significantly more 
sensitive to approximations than other quantities such 
as the internal energy. 



1. Pair Density Matrix 

If one only used the bare potential as in Eq. ^ (the 
primitive approximation for the action) , the convergence 
would be very slow pij and would result in an extremely 
inefficient many-particle simulations. Instead, we nu- 
merically solve the two-particle problem with the matrix 
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FIG. 2: Accuracy study of PIMC simulations of the isolated 
hydrogen atom using different orders n\ in the expansion for- 
mula Fl2l for the action and energy. The calculated 2K + V 
(exact value equals zero) and the deviation from the exact 
potential energy of —27.2 eV are shown for different orders 
from simulations at T = 10 000 K using r" 1 = 10 6 K. 



squaring technique |15j . Numerical representations of the 
exact pair density matrices are stored in tables used by 
the PIMC simulation program by expanding in the small 
variables s and z: 

u(r,r';r) = - [u (r; r) + u (r'; r)] 



HA k 

^^ Ufe ,(g;r)z 2 ^ 2 ^), (12) 
fc=i 3=0 



where 



|r'|) a = 



(13) 



and r and r' denote the separation of the two particles 
at adjacent time slices. The accuracy of these tables is 
crucial for all computed results. Using the precomputed 
pair density matrices allows one to employ a much larger 
time step because one starts with a solution of the two- 
particle problem. Fig. || shows how accurate this method 
is. The internal energy of an isolated hydrogen atom at 
sufficiently low temperature (T = 10 000 K) in a large box 
(L = 26) is compared with the exact ground state energy 
of — 13.6eV. The temperature was chosen low enough so 
that excited states can be neglected; the contribution to 
the energy from the occupation of the first excited state 
is 7- 10 -5 eV at this temperature. Also shown is how well 
the kinetic energy K and the potential energy V satisfy 
the virial theorem 2K + V = 0, thus determining the ac- 
curacy with which the pressure can be determined. For 
a time step of r _1 = 10 6 K, the analysis shows a quick 
convergence with the order of terms considered in the ac- 
tion expansion Eq. (h2). Using terms up ua — 3 reduces 
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the error to 0.033 (3) eV in energy and to 0.039 (8) eV 
for 2K + V. This corresponds to an inaccuracy in the 
pressure equivalent to a non-interacting molecular gas at 
T = 260 (30) K. 



2. Time step dependence 




t" 1 (10 6 K) 

FIG. 3: The error in the virial as a function of the number of 
time slices for an isolated hydrogen molecule. The classical 
nuclei are hxed at the equilibrium bond length, thus the exact 
value for 2K + V is zero at a sufficiently low temperature. 



Employing the pair approximation of the density ma- 
trix does not include correlation effects for three or more 
particles (for example between one electron and two pro- 
tons). We now estimate how small the time step must 
be to obtain a given accuracy for an isolated molecule. 
Fig. U shows results for different time steps and tem- 
peratures with the nuclei kept fixed at the equilibrium 
position of R = 1.4008. From the virial theorem, it fol- 
lows that 2K + V = at a sufficiently low temperature. 
The exact energy per atom is — 15.973eV [ pH - The T 
dependence is small suggesting that the electrons are in 
the ground state. However, one finds a significant depen- 
dence on the time step. Using r 1 = 2 x 10 6 K reduces 
the error in the energy per atom to 0.036(3) eV and in 
2K + V to 0.090 (16) eV. The time step error is larger 
than the errors of the inaccuracies in the pair density 
matrices discussed above. The error in 2K + V corre- 
sponds to a pressure of an non-interacting molecular gas 
at T = 700 (100) K, which provides us with an approxi- 
mate limit of accuracy in the equation of state calculated 
from many-particle simulations with r _1 = 2 x 10 6 K. 



3. Finite Size Dependence 

The estimation of the finite size errors is more diffi- 
cult to assess because the needed PIMC simulations are 
computationally much more demanding. The required 
computer time increases rapidly with the number of par- 
ticles making it challenging to obtain converged results 
for paths corresponding to large systems. 

Most results from many-body simulations reported in 
this work were calculated with N — 32 pairs of elec- 
trons and protons in a periodically repeated simulation 
cell. To study the effect on N, we performed simulations 
for N — 16 and 64 pairs of protons and electrons for a 
density of r s = 2.6 and T > 10 000 K. We chose the 
highest density under consideration because one expects 
the finite size dependence to be largest there due to the 
stronger interaction between the atoms. 
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FIG. 4: Finite size error of the pressure as a function of tem- 
perature relative to simulations with N = 64: pairs of protons 
and electrons at a density of r s = 2.6. 



The finite size dependence of the pressure, shown in 
Fig. |J, is small at high temperatures but grows to ap- 
proximately 4 (2) % near T — 30 000 K. In this regime, 
the hydrogen undergoes structural changes involving the 
formation of atoms, which affect the pressure. This study 
provides us only with an estimate of the finite size de- 
pendence. An extrapolation to N — * oo would require 
significantly larger systems, not currently feasible at low 
temperatures. 

Fig. U shows the finite size error of the internal en- 
ergy The smaller systems are more strongly bound by 
approximately 0.2 eV per atom, probably because of the 
interaction of a charge with its own image. For lower 
densities, we expect this value to be smaller. 
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FIG. 5: Finite size error of the internal energy as a function 
of temperature relative to simulations with ./V = 64 pairs of 
particles at a density of r 3 — 2.6. 




FIG. 6: The effect on the pressure of two different nodal sur- 
faces: of the free particle density matrix and of the variational 
density matrix. 



4- Nodal Approximation 

In the above, we have studied controlled approxima- 
tions. The only uncontrolled approximation in the re- 
stricted PIMC method is the use of trial density matrix 
to constrain the paths. The nodal surfaces are impor- 
tant only if the electrons are degenerate: at low temper- 
atures or at high densities. Recall that in this work we 
focus on hydrogen only at low density, where the elec- 
trons are bound in atoms and molecules and have a low 
or moderate degeneracy. Even at low density, one still 
needs a nodal surface in order to prevent the formation 
of unphysical clusters like H$ and H4 or even the col- 
lapse of the entire system, but the precise shape of the 
nodes is not important at low density as shown in Fig. |^. 
Comparing FP and VDM nodes for r s = 2.6, one only 
finds differences in the pressure for T < 15 625 K, which is 
approximately where, at this density, the syst em sh ows 
a significant molecular signature (see section IIIB and 
Fig. |l6|). In this regime, FP nodes systematically lead to 
a too high pressure, while simulations with VDM nodes 
stay closer to the prediction of a semi-empirical chemical 
model At a lower density, r s = 4, results from FP 
and VDM nodes agree within the error bars. The dif- 
ferences in the internal energy, as shown in Fig. [?], using 
FP and VDM nodes are significantly smaller than the 
pressure deviations. One either finds agreement within 
the error bars or that VDM nodes predict lower inter- 
nal energies, which was used in |L7| to show that VDM 
nodes are the more accurate nodal surface. The observed 
energy differences did not exceed 0.1 eV per atom. 

For even lower density, the nodes are less relevant be- 
cause they become only important in a collision of two 
molecules, which occur less frequently at lower density. 
This trend can also be understood in terms of the degen- 
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FIG. 7: Internal energy computed with PIMC using two dif- 
ferent nodal surfaces: the free particle density matrix and the 
variational density matrix. 



eracy of the electrons. The degree of degeneracy mani- 
fests itself in the path integral formalism by the probabil- 
ity for the electrons to be involved in a permutation. At 
high temperature, the paths are very short and permu- 
tations are rare. At low temperature and high density, 
the paths are long and can form long permutation cy- 
cles. However in hydrogen at low density, the paths are 
localized due to the attraction in atoms and molecules 
and permutations are rare. Fig. [I] shows that the per- 
mutation probability never reaches 1% for r s = 4 (see 
Tab. |lj). For higher densities, the permutation probabil- 
ity is increased as indicated by the contour lines. This 
is consistent with the temperature and density depen- 
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dence of the nodal error in pressure and internal energy 
discussed above. 



III. RESULTS 



A. Equation of State 
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FIG. 8: Internal energy per atom vs. temperature for a den- 
sity of r s = 10 comparing the SC-EOS jj] with PIMC calcu- 
lations. 
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FIG. 9: Internal energy per atom vs. temperature for a den- 
sity of r s — 10 as shown in Fig. |^ but here for lower temper- 
atures also including results from the activity expansion J18[ . 

Tab. | gives the complete set of energies and pressures 
at 6 densities and 8 temperatures. We now compare these 
results with several models for hydrogen. We begin our 
discussion by studying the internal energy per atom as 
a function of temperature shown in Figs. |^ - [ll] for two 
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FIG. 10: Internal energy per atom vs. temperature as shown 
in Fig. H but here for a density of r s = 2.6. 
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FIG. 11: Internal energy per atom vs. temperature for a den- 
sity of r 3 = 2.6 as shown in Fig. |l^ but here for lower tem- 
peratures also including results from the activity expansion 
(ActEx) (li) and the fluid variational theory (FVT) M. 



selected densities corresponding to r s — 10.0 and 2.6. 
Generally, we find a fairly good overall agreement with 
the EOS by SC Q over the entire temperature and den- 
sity range discussed in this work. The agreement is par- 
ticularly good in the molecular and atomic regime for 
r a = 10.0, as shown in Fig. |. There the SC energies are 
within the error bars of the PIMC results. At higher tem- 
perature shown in Fig. |^, we find systematic deviations 
of up to 5 eV per atom at T = 250 000 K. They indicate 
that the SC energies are too low at high temperatures 
and too high at intermediate temperatures (see Fig. ||). 
One possible explanation for the deviations at high tem- 
perature is that the SC model underestimates the degree 
of ionization (see discussion in pfj ) . 
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We also studied these deviations as a function of den- 
sity. The cross-over temperature, above which the SC- 
EOS underestimates the energy, increases with density. 
At r s — 10.0, the cross-over is near 70 000 K compared to 
130 000 K (Fig. at r s = 2.6. At temperatures below 
20 000 K for r s = 2.6, one also finds some small deviations 



up to 0.5 eV per atom (Fig. 11) 
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FIG. 13: Pressure vs. temperature at a density of r a = 10 
showing results from the fluid variational theory [j^], the SC- 
EOS J5J, and PIMC simulations. 



tures 62 500 > T > 15 625 K, one finds pressure differ- 
ences, which are of the same magnitude as the finite size 
effects in PIMC. For temperatures below 15 625 K, the 
increased statistical errors in the PIMC pressure are of 
the same size as the observed deviations. 



FIG. 12: Internal energy per atom vs. density for different 
temperatures from SC-EOS (jj, the activity expansion (Ac- 
tEx) jn| (not shown for 5000 K since nearly identical to SC) 
and PIMC calculations . 

Fig. [12] shows a comparison of energy vs. density for 
several temperatures. It shows that the SC-EOS over- 
estimates the energy for T = 5000 K and 31 250 K and 
underestimates it for 125 000 K for densities higher than 
those corresponding to r s = 2.6. 

Now let us compare the pressure from the SC-EOS 
with that from the PIMC simulation using Eq. |ll] in 
Tab. Q. We find remarkably good agreement of the entire 
range of temperature and density under consideration. 
For a low density such as r s = 10.0, this is shown in 
Fig. [ll| As expected, one finds that both methods in- 
terpolate between the limit of an ideal Fermi gas at high 
temperatures and non-interacting molecular gas at low 
temperatures. Fig. [l4| confirms the good agreement at a 
higher density of r s = 2.6. As a result of the strong inter- 
actions at this density, one finds that the pressure at low 
temperatures is significantly above the non-interacting 
molecular gas limit. 

We find that the SC-EOS underestimates the pressure 
by about 3% for T > 62 500 K. This difference is outside 
the error b ar fr om the approximations in PIMC discussed 
in section [I B and could be interpreted as a further in- 
dication, in addition to the observed energy deviations, 
that the SC model underestimates the degree of ioniza- 
tion at high temperatures. For intermediate tempera- 
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FIG. 14: Pressure vs. temperature at a density of r s 
shown in Fig. h_3l 
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In Fig. EsL we show the pressure as a function of den- 
sity, which confirms the good agreement. The figure also 
indicates that, at 5000 K and r s > 4, the pressure is close 
to the pressure of a non-interacting molecular gas. 

In our comparison, we also included results from the 
activity expansion by Rogers ]l8fl , which shows very good 
agreement in pressure and internal energy (see Figs. |{], 
|ll], and ||). The differences are small but increase with 
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FIG. 15: Pressure vs. density for different temperatures. 



density. In the molecular and atomic regime, one also 
finds good agreement with the fluid variational theory 
(FVT) by Juranek and Redmer as shown in Figs. |], 
[y], and [l4|. For higher temperatures, the FVT model 
is not applicable since it does not include ionization of 
atoms. 



B. Pair Correlation Functions 

There are four different pair correlation functions 
which can be directly obtained from many-body simu- 
lations and provide direct information about the state of 
the system. Shown in the following figures is an extensive 
set of pair correlations which allow one to estimate the 
microscopic structure of the system and allow a direct 
comparison with other simulations. The proton-proton 
pair correlation functions from PIMC simulations with 
free particle nodes are shown in Fig. [16[ For T ~ 20 000 K 
a peak at the bond length of 1.4 emerges, which clearly 
demonstrates the formation of molecules. We found it 
useful to multiply the pair correlation function by an ex- 
tra density factor n = N/v so that the area under the 
peak is proportional to the molecular fraction. The peak 
height gets smaller with decreasing density as a result of 
entropic dissociation of the molecules, driven by the num- 
ber of unbound states at low density. Thermal dissocia- 
tion also reduces the number of molecules with increasing 
temperature. For r s ~ 2, we expect that pressure dissoci- 
ation diminishes the number of molecules with increasing 
density jl4| but this density range is beyond the scope of 
this paper. 

The proton-electron pair correlation function multi- 
plied by the density is shown in Fig. [Tt]. The peak 
near the origin shows the increased probability of finding 
an electron near a proton due to the Coulomb attrac- 
tion. The peak height decreases with temperature and 



TABLE I: Pressure and internal energy per atom and the 
resulting Hugoniot from PIMC simulations with 32 pairs of 
particles and r _1 = 2 • 10 6 K using free particle nodes except 
for * where VDM nodes were employed instead. The prob- 
abilities x for finding a proton in a given state for the three 
dominant species are derived from a cluster analysis. P pe rm 
is the permutation probability for the electrons. 



T(K) P(GPa) 



E(eV) x H+ xu xn 2 P P 



14.0 250 000 



14.0 
14.0 
14.0 
14.0 
14.0 
14.0 
14.0 



125 000 
62 500 
31250 
15 625 
10 000 
7812 
5 000 



4.000 (2 
1.955 (2 
0.901 (2 
0.334 (2 
0.127(2 
0.081 (4 
0.047 (5 
0.028 (6 



10.0 250 000 
10.0 125 000 
62 500 



10.0 
10.0 
10.0 
10.0 
10.0 
10.0 



31250 
15 625 
10 000 
7812 
5 000 



10.902 (4 
5.259 (6 
2.329 (5 
0.831 (5 
0.344 (4 
0.198 (9 
0.144(7 
0.068(15) 



6.0 250 000 



6.0 
6.0 
6.0 
6.0 
6.0 
6.0 
6.0 



125 000 
62 500 
31250 
15 625 
10 000 
7812 
5 000 



49.46 (3 
23.00 (3 
9.56 (2 
3.58 (3 
1.52(2 
0.77(4 
0.63 (5 
0.29 (9 



4.0 250 000 
4.0 125 000 
62 500 
31250 
15 625 
10 000 
7812 
5 000 



4.0 
4.0 
4.0 
4.0 
4.0 
4.0 



162.46(10 
73.00 (23 
29.75(16 
11.22(22 
5.01(17 
3.23 (30 
2.20(14 
1.19 (25 



3.0 250 000 
3.0 125 000 
62 500 
31250 
15 625 
10 000 



3.0 
3.0 
3.0 
3.0 



374.47(14 
165.21 (22 
67.70 (25 
26.67 (24 
13.08 (28 
9.26 (26 



2.6 250 000 
2.6 125 000 
2.6 62 500 
2.6 31250 
2.6 15 625 
2.6 10 000 



566.4 (4 
246.0 (5 
101.7(4 
43.1 (5 
19.0 (8 
11.7(9 



62.93 (4) 
29.97(4) 
11.85 (4) 
-2.97(3) 
-11.30 (4) 
-12.43 (6) 
-13.34(13) 
-15.00(12) 



0.99 0.01 0.00 

0.95 0.05 0.00 

0.64 0.36 0.00 

0.24 0.76 0.00 

0.12 0.72 0.15 

0.08 0.58 0.33 

0.01 0.11 0.88 



0.000 
0.000 
0.000 
0.000 
0.000 
0.000 
0.000 
0.000 



62.00 (3) 
28.56 (3) 
9.41 (3) 
-4.91 (3) 
-11.49 (3) 
-12.79(5) 
-13.54 (9) 
-15.05 (7) 



0.98 0.02 0.00 

0.90 0.10 0.00 

0.57 0.43 0.00 

0.23 0.74 0.02 

0.11 0.60 0.27 

0.03 0.34 0.62 

0.00 0.07 0.92 



0.000 
0.000 
0.000 
0.000 
0.000 
0.000 
0.000 
0.000 



59.33 (4) 
24.70 (4) 
4.79 (3) 
-6.92 (4) 
-11.83 (4) 
-13.38 (6) 
-13.95 (7) 
-15.17(12) 



0.94 0.06 0.00 

0.80 0.20 0.00 

0.52 0.46 0.01 

0.21 0.68 0.09 

0.08 0.49 0.40 

0.04 0.38 0.56 

0.00 0.09 0.90 



0.000 
0.000 
0.000 
0.000 
0.001 
0.000 
0.000 
0.000 



55.63 (4) 
20.24 (9)* 
1.23 (6)* 
-8.32 (8)* 
-11.87 (6)* 
-13.43(11)' 
-14.29 (6) 
-15.20 (9) 



0.88 0.12 0.00 

0.72 0.26 0.00 

0.47 0.46 0.03 

0.19 0.58 0.18 

0.03 0.30 0.63 

0.01 0.18 0.80 

0.00 0.11 0.88 



0.000 
0.000 
0.001 
0.004 
0.008 
0.005 
0.004 
0.002 



51.79(2) 
16.58 (4) 
-1.04 (4) 
-9.02 (4) 
-12.29 (4) 
-13.79 (4) 



0.83 0.15 0.01 

0.66 0.29 0.01 

0.45 0.43 0.06 

0.15 0.42 0.34 

0.03 0.35 0.60 



0.000 
0.001 
0.005 
0.028 
0.059 
0.037 



49.58 (4) 
14.57 (5)* 
-2.25 (4)* 
-9.38 (5)* 
-12.51 (9)* 
-13.73(10)* 



0.80 0.17 0.01 
0.65 0.28 0.02 
0.42 0.38 0.09 
0.14 0.42 0.33 
0.02 0.34 0.60 



0.000 
0.002 
0.014 
0.078 
0.168 
0.126 



increases with density because of thermal ionization and 
entropy ionization respectively. At low temperature, the 
peak can be interpreted as occupation of bound states al- 
though (unbound) scattering states can also contribute. 
From proton-electron pair correlation alone, one cannot 
distinguish between an atomic and a molecular state. 

Figs. [l8] and |l9| show the electron-electron pair cor- 
relation functions for pairs with anti-parallel spin. The 
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FIG. 16: Proton-proton pair correlation function multiplied 
by the density n. The columns correspond to to r s values and 
the rows to different temperatures T. 



FIG. 18: Electron-electron pair correlation function for elec- 
trons with opposite spin multiplied by the density n is shown. 
The columns correspond to different r s values and the rows 
to different temperatures T. 



T(K) 



62500 



31250 



15625 



r =14 r =10 r=6 r=4 r=3 r=2.6 



n g„ e (r) 



10000 



7812 



5000 



0240240240240240246 

r 

FIG. 17: Proton-electron pair correlation function multiplied 
by the density n. The columns correspond to different r s 
values and the rows to different temperatures T. 



peak at small separations comes from the formation of 
the molecular bond. For pairs of electrons with parallel 
spin, one always finds a strong repulsion due to the Pauli 
exclusion principle and to a lesser extent to the Coulomb 
repulsion. This is shown in Fig. 19. 



C. State of Hydrogen 

In this section, we discuss the phase diagram of hy- 
drogen as shown in Fig. |l|. The diagram shows the ap- 
proximate location of the molecular, the atomic and the 
plasma regimes. The PIMC simulations, since they are 
based on the basic description in terms of electrons and 
protons, do not directly lend themselves to determining 
the number of compound particles such as molecules and 
atoms. (Methods for determining this from PIMC sim- 
ulations will be discussed in a future publication.) In 
order to obtain an estimate of the atomic and molecu- 
lar fraction, we employed a cluster anal ysi s of the PIMC 
path configurations. As described in |2lf| , we consider 
two protons as belonging to one cluster if they are less 
than 1.9 ae apart. An electron belongs to one particular 
cluster if it is less than 1.4 ob away from any proton in 
the cluster. The two cut-off radii were chosen from the 
molecular and atomic ground state distribution. This 
analysis gives reasonable estimates for the molecular and 
atomic fractions at low temperatures. At high tempera- 
ture, it overestimates the number of bound states because 
even in an (unbound) collision, two particles are counted 
falsely as being part of the cluster. We corrected for this 
by applying an additional criterion: a particle can only 
be considered as bound if the difference in action to re- 
move it from the cluster is positive. This method leads 
to the expected corrections at high temperature. (The 
regime boundaries in Fig. [j] discussed below are hardly 
affected by the additional correction.) 

The lower dashed line represents the region where 60% 
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T(K) 



62500 



r =14 r =10 r=6 r=4 r=3 r=2.6 




10 15 



FIG. 19: Electron-electron pair correlation function. The 
solid lines correspond to pairs of electrons with parallel spin. 
For the sake of comparison, we also show the pair correlation 
functions of pairs with opposite spin as dashed lines. This 
function is strongly peaked near the origin in the presence of 
molecules as shown in Fig. [l^. The columns correspond to 
different r s values and the rows to different temperatures T. 



of the protons are bound in molecules. When the number 
of protons bound in atoms ( i. e. with an electron) drop 
below 40% we labeled this state as plasma as shown in the 
upper dashed-dotted line. It should be emphasized that 
the location of these lines depends on the choice of these 
limits as well as on the cut-off radii used to determine 
the clusters in this place. Fig. [j] also shows the location 
of isobars, which appear as almost straight lines in this 
double logarithmic graph. The slope is different from the 
ideal gas because the pressure depends on ionization and 
dissociation. 

Tab. I shows the fraction of the three most frequently 
found species: molecules, atoms and free protons where 
x is defined as the probability of finding a proton in a 
certain compound particle. It should be noted that the 
sum: a; H +-|-a;H+^H2 is l ess than 1 since other clusters 
have a non-zero probability. The largest contributions 
besides those listed are with a maximum of 0.06 for 
r s = 2.6 and T = 15 625 K followed by H 3 with x < 0.03 
and H~ with x < 0.02. Even larger clusters occur very 
infrequently. The cluster analysis also gives an estimate 
for the fraction of free electrons, which agrees well with 
the number of ionized protons: ■ 

Fig. po| shows a comparison of the fraction of molecules 
and ionized atoms for r s — 10. One finds that the molec- 
ular fractions decays rapidly with temperature. The re- 
sulting atoms are then ionized at even higher temper- 
atures leading to the observed increase in the number 
of free protons. The PIMC predictions for the molecu- 
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FIG. 20: Fraction of molecules and free protons as a function 
of temperature for r s = 10 comparing the cluster analysis of 
the PIMC results with the SC model Q and the fluid varia- 
tional theory (FVT) E§. 
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FIG. 21: Fraction of molecules and free protons as shown in 
Fig. ^ but here for a density of r s = 2.6. 



lar fraction agree very well with the SC model as well 
as with the FVT. On the other hand, the PIMC results 
shows a significantly higher degree of ionization than SC. 
The same comparison for a higher density of r s = 2.6 in 
Fig. ||l| shows that the cluster analysis leads to smaller 
number of molecules than predicted by the SC model. 

Summarizing one can say that the cluster analysis pro- 
vides us with reasonable estimates for the number of 
atoms and molecules in the considered density range. We 
caution the reader that other definitions of atoms and 
molecules, possible in a many-body hydrogen, while giv- 
ing qualitatively similar results, may differ quantitatively. 
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Our computed numerical values must be used with cau- 
tion. A rigorous, many-body definition of a bound state 
remains to be applied. Several ideas are discussed in p2| . 

IV. CONCLUSIONS 

In this work, we studied the high-temperature equation 
of state of hydrogen at low and intermediate densities 
and find a remarkably good agreement with the SC-EOS. 
Generally, one finds that the deviations in the energy are 
more pronounced than the differences in the pressure. 
We find significant deviations in the EOS of temperatures 
« lOOOOOi-T, most likely caused by an underestimate of 
the degree of ionization at those temperatures. In future 
work, we will extend this comparison to higher densities. 
There one expects to find substantial differences between 

I 



the SC model and PIMC, which manifest themselves in 
a different shock Hugoniot JT^j . 
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